PROJEKT 2 Predykcja zanieczyszczenia PMx

KRZYSZTOF KOWALSKI 407142

Zadanie polega na przeprowadzeniu analizy i predykcji szeregów czasowych dotyczących poziomu pyłu zawieszonego przy użyciu modelu XGBoost. Pył PM10 i PM2.5 to drobne cząstki stałe lub ciekłe zawieszone w powietrzu, których średnica wynosi odpowiednio 10 lub mniej i 2,5 mikrometra lub mniej.

Zanieczyszczenie powietrza tymi cząsteczkami stanowi istotny problem dla zdrowia publicznego, ponieważ może przenikać głęboko do płuc i układu krwionośnego, co ma negatywny wpływ na zdrowie. Predykcja zawartości PMx w powietrzu pozwoli na lepsze planowanie i informowanie ludności o potencjalnych zagrożeniu.

Część 1: Przygotowanie danych i inżynieria cech (20 punktów)

1. Tworzenie cech (15 punktów)

Utwórz cechy (features) w oknach czasowych oraz cechy oparte na opóźnieniach wartości pyłu PM2.5 w celu uwzględnienia historii w predykcji.

Dodatkowo, wykorzystaj informacje z datatime do tworzenia dodatkowych cech, takich jak dzień tygodnia, miesiąc, rok itp

Uzasadnij swój wybór.

Czy jest możliwe uwzględnienienie czynnika przestrzennego w tworzeniu cech?

Wczytanie danych

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from statsmodels.tsa.seasonal import STL
import xgboost as xgb
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import mean_squared_error, r2_score
from statsmodels.graphics.tsaplots import plot_acf
import statsmodels.api as sm
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import TimeSeriesSplit

Po zapoznaniu sie z wygladem danych w excelu zdecydowano na wcztanie sie ich w taki sposob aby kolumnami byly stacje pomiarowe natomiast indeksem data pomiaru.

data = pd.read_excel('2019_PM25_1g.xlsx', skiprows = 5)
data = data.melt(id_vars=["Kod stanowiska"], var_name="station_code", value_name="PM25").rename(columns = {'Kod stanowiska':'Data', 'station_code':'Nazwa Stacji'}) 
data['Nazwa Stacji'] = data['Nazwa Stacji'].str[:-9]
data

	Data	Nazwa Stacji	PM25
0	2019-01-01 01:00:00	DsDusznikMOB	33.40530
1	2019-01-01 02:00:00	DsDusznikMOB	13.80280
2	2019-01-01 03:00:00	DsDusznikMOB	9.94056
3	2019-01-01 04:00:00	DsDusznikMOB	6.75889
4	2019-01-01 05:00:00	DsDusznikMOB	7.88722
...	...	...	...
551875	2019-12-31 20:00:00	ZpSzczPils02	19.16890
551876	2019-12-31 21:00:00	ZpSzczPils02	19.71360
551877	2019-12-31 22:00:00	ZpSzczPils02	20.24390
551878	2019-12-31 23:00:00	ZpSzczPils02	24.74170
551879	2020-01-01 00:00:00	ZpSzczPils02	36.85750

551880 rows × 3 columns

W celu lepszej ekslopracji i poznania zbioru stwierdzono, ze nalezy dane uzupelnic o lokalizacje stacji. Informacje te są łatwo dostępne pod adresem https://powietrze.gios.gov.pl/pjp/archives/.

stacje = pd.read_excel('kody_stacji.xlsx')
stacje

	Nr	Kod stacji	Kod międzynarodowy	Nazwa stacji	Stary Kod stacji \n(o ile inny od aktualnego)	Data uruchomienia	Data zamknięcia	Typ stacji	Typ obszaru	Rodzaj stacji	Województwo	Miejscowość	Adres	WGS84 φ N	WGS84 λ E
0	1	DsBialka	NaN	Białka	NaN	1990-01-03	2005-12-31	przemysłowa	podmiejski	kontenerowa stacjonarna	DOLNOŚLĄSKIE	Białka	NaN	51.197783	16.117390
1	2	DsBielGrot	NaN	Bielawa - ul. Grota Roweckiego	NaN	1994-01-02	2003-12-31	tło	miejski	w budynku	DOLNOŚLĄSKIE	Bielawa	ul. Grota Roweckiego 6	50.682510	16.617348
2	3	DsBogatFrancMOB	PL0602A	Bogatynia Mobil	DsBogatMob	2015-01-01	2015-12-31	tło	miejski	mobilna	DOLNOŚLĄSKIE	Bogatynia	ul. Francuska/Kręta	50.940998	14.916790
3	4	DsBogChop	PL0315A	Bogatynia - Chopina	NaN	1996-01-01	2013-12-31	przemysłowa	miejski	kontenerowa stacjonarna	DOLNOŚLĄSKIE	Bogatynia	ul. Chopina 35	50.905856	14.967175
4	5	DsBogZatonieMob	PL0576A	Bogatynia - Mobil	NaN	2012-01-01	2012-12-31	przemysłowa	miejski	mobilna	DOLNOŚLĄSKIE	Bogatynia	ul. Konrada, Zatonie	50.943245	14.913327
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
1090	1091	ZpSzczPilsud	PL0249A	Szczecin, ul. Piłsudskiego	ZpSzczecin002, ZpSzczPils02	2004-12-31	NaT	komunikacyjna	miejski	kontenerowa stacjonarna	ZACHODNIOPOMORSKIE	Szczecin	ul. Piłsudskiego 1	53.432169	14.553900
1091	1092	ZpSzczWSSEEnerg	NaN	Energetyków	NaN	1992-01-01	2003-12-31	komunikacyjna	miejski	w budynku	ZACHODNIOPOMORSKIE	Szczecin	ul. Energetyków 2	53.420475	14.561934
1092	1093	ZpSzczWSSESped6	NaN	Spedytorska	NaN	1992-01-01	2004-01-01	tło	miejski	kontenerowa stacjonarna	ZACHODNIOPOMORSKIE	Szczecin	ul. Spedytorska 6	53.415043	14.555347
1093	1094	ZpWalWalczWSSE	NaN	Wałcz	NaN	1992-01-01	2005-01-01	tło	miejski	w budynku	ZACHODNIOPOMORSKIE	Wałcz	ul. Bydgoska 86	53.263667	16.492596
1094	1095	ZpWiduBulRyb	PL0182A	Widuchowa	ZpGryfWiduchowa003	2003-01-01	NaT	tło	pozamiejski	kontenerowa stacjonarna	ZACHODNIOPOMORSKIE	Widuchowa	ul. Bulwary Rybackie 1	53.122325	14.382245

1095 rows × 15 columns

Ponieważ niektóre z naz stacji zmieniły nazwę, należało je odpowiednio zamienić.

data.rename(columns = {"LbNaleczow": "LbNaleczAlMa", "MzKonJezMos": "MzKonJezWieMOB", "PdSuwPulaskp": "PdSuwPulask2", "PmGdaLeczk08": "PmGdaLeczkow", "ZpSzczAndr01": "ZpSzczAndrze", "ZpSzczPils02": "ZpSzczPilsud"}, inplace = True)
data = data.merge(stacje, left_on = 'Nazwa Stacji', right_on = 'Kod stacji', how = 'inner')
data = data[['Data', 'Nazwa stacji', 'PM25', 'Typ obszaru', 'Województwo', 'WGS84 φ N', 'WGS84 λ E']]
data

	Data	Nazwa stacji	PM25	Typ obszaru	Województwo	WGS84 φ N	WGS84 λ E
0	2019-01-01 01:00:00	Duszniki-Zdrój	33.40530	miejski	DOLNOŚLĄSKIE	50.402645	16.393319
1	2019-01-01 02:00:00	Duszniki-Zdrój	13.80280	miejski	DOLNOŚLĄSKIE	50.402645	16.393319
2	2019-01-01 03:00:00	Duszniki-Zdrój	9.94056	miejski	DOLNOŚLĄSKIE	50.402645	16.393319
3	2019-01-01 04:00:00	Duszniki-Zdrój	6.75889	miejski	DOLNOŚLĄSKIE	50.402645	16.393319
4	2019-01-01 05:00:00	Duszniki-Zdrój	7.88722	miejski	DOLNOŚLĄSKIE	50.402645	16.393319
...	...	...	...	...	...	...	...
499315	2019-12-31 20:00:00	Szczecin_Budzysza Wosia	15.59590	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294
499316	2019-12-31 21:00:00	Szczecin_Budzysza Wosia	13.18250	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294
499317	2019-12-31 22:00:00	Szczecin_Budzysza Wosia	14.37640	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294
499318	2019-12-31 23:00:00	Szczecin_Budzysza Wosia	17.79410	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294
499319	2020-01-01 00:00:00	Szczecin_Budzysza Wosia	20.86120	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294

499320 rows × 7 columns

EDA

Zaczęto od sprawdzenia ilości wartości brakujących w każdej z kolumn. Widoczne są wartości brakujące w kolumnie związanej z wielkością zanieczyszczenia. W celu lepszej eksploracji tych wartości sprawdzono jak prezentują się dane w tej kolumnie w zależności od stacji z której pochodzi informacja.

for col in data.columns:
    print(f'W kolumnie {col} znajduje się {np.sum(data[col].isna())} co stanowi {round(100 * np.sum(data[col].isna()) / len(data[col]), 2)}% wszystkich wartości.')

W kolumnie Data znajduje się 0 co stanowi 0.0% wszystkich wartości.
W kolumnie Nazwa stacji znajduje się 0 co stanowi 0.0% wszystkich wartości.
W kolumnie PM25 znajduje się 18660 co stanowi 3.74% wszystkich wartości.
W kolumnie Typ obszaru znajduje się 0 co stanowi 0.0% wszystkich wartości.
W kolumnie Województwo znajduje się 0 co stanowi 0.0% wszystkich wartości.
W kolumnie WGS84 φ N znajduje się 0 co stanowi 0.0% wszystkich wartości.
W kolumnie WGS84 λ E znajduje się 0 co stanowi 0.0% wszystkich wartości.

fig, axs = plt.subplots(8, 8, figsize=(16, 16))

for i, col in enumerate(data['Nazwa stacji'].unique()):
    ax = axs[i // 8, i % 8]
    ax.hist(data[data['Nazwa stacji'] == col]['PM25'], bins = 20, color = 'skyblue', edgecolor = 'black')
    ax.set_title(col)
    ax.set_xlabel('Wartość')
    ax.set_ylabel('Częstotliwość')
    ax.grid()

plt.tight_layout()

Stacje posiadają podobny rozkład częstotliwości wystąpień. Jest to rozkład nie przypominający rozkładu normalnego.

brakujace_wartosci = data[data['PM25'].isnull()].groupby('Nazwa stacji')['PM25'].size()
brakujace_wartosci

Nazwa stacji
Białystok, ul. Waszyngtona                217
Bielsko-Biała, ul.Partyzantów              20
Borsukowizna, Szkółka Leśna                44
Bydgoszcz, pl. Poznański                  849
Duszniki-Zdrój                            114
Dębica, ul.Grottgera                      311
Elbląg, ul. Bażyńskiego                  1055
Horyniec-Zdrój-Park                       673
Jarosław, ul. Pruchnicka                   43
Jawor                                      49
Jelenia Góra, ul. Ogińskiego              262
Kalisz, ul. Wyszyńskiego                  212
Katowice, ul. Kossutha                    569
Kielce, ul. Targowa                       419
Kraków, Aleja Krasińskiego                 30
Kraków, ul. Bujaka                        100
Kraków, ul. Bulwarowa                      77
Kędzierzyn-Koźle, ul. Śmiałego            590
Legionowo, ul. Zegrzyńska                 123
Lublin, ul. Obywatelska                   348
Mielec, ul. Biernackiego                  383
Mińsk Mazowiecki-Kazikowskiego            137
Mogilno,Kościuszki                        441
Nowa Sól, ul. T. Kościuszki               596
Olsztyn, ul. Puszkina                       2
Otwock, ul. Brzozowa                       57
Piastów, ul. Pułaskiego                   228
Prudnik mobilna                           128
Przemyśl, ul. Grunwaldzka                 515
Płock, ul. Reja                            77
Radom, ul. Tochtermana                    151
Rymanów Zdrój, ul. Parkowa                224
Rzeszów, ul. Piłsudskiego                  11
Siedlce, ul. Konarskiego                  435
Skarżysko-Kamienna, MOBILNA               801
Starachowice, ul. Złota                    41
Szczecin_Budzysza Wosia                  1620
Tarnów, ul. Ks. Romana Sitko               12
Toruń, ul. Dziewulskiego                 1547
WIOŚ Gołdap-Uzdrowisko                    607
Warszawa, al. Niepodległości              205
Warszawa, ul. Bajkowa                      26
Warszawa, ul. Chrościckiego                49
Warszawa, ul. Kondratowicza                51
Warszawa, ul. Tołstoja                    454
Warszawa, ul. Wokalna                     551
Wrocław, al. Wiśniowa                      53
Wrocław, wyb. Conrada-Korzeniowskiego     114
Wschowa ul. Kazimierza Wielkiego          497
Włocławek, ul. Okrzei                      47
Zgierz, ul. Mielczarskiego                492
Zielona Góra ul. Krótka                   622
Złoty Potok, Leśniczówka                   94
Łódź, ul. Czernika                        114
Łódź, ul. Gdańska                         677
Żary, ul. Szymanowskiego                  264
Żyrardów, ul. Roosevelta                  232
Name: PM25, dtype: int64

Wyróżniają się stacje posiadające większą ilość wartości brakujących niż inne. Może to wynikać z innych częstotliwości próbkowania niż co godzinę, chwilowej awarii lub wyłączenia stacji lub pózniejszym pojawieniu się stacji niż badany okres.

Stwierdzono, że odpowiednim uzupełnieniem wartości brakujących będzie wybranie średnie z tej samej daty dla stacji znajdujących się w tym samym województwie.

data['PM25'] = data['PM25'].fillna(data.groupby(['Data', 'Województwo'])['PM25'].transform('mean'))
brakujace_wartosci = data[data['PM25'].isnull()].groupby('Nazwa stacji')['PM25'].size()
brakujace_wartosci

Nazwa stacji
Białystok, ul. Waszyngtona           7
Borsukowizna, Szkółka Leśna          7
Bydgoszcz, pl. Poznański             1
Kalisz, ul. Wyszyńskiego           212
Kielce, ul. Targowa                  1
Kędzierzyn-Koźle, ul. Śmiałego       1
Legionowo, ul. Zegrzyńska            1
Lublin, ul. Obywatelska            348
Mińsk Mazowiecki-Kazikowskiego       1
Mogilno,Kościuszki                   1
Otwock, ul. Brzozowa                 1
Piastów, ul. Pułaskiego              1
Prudnik mobilna                      1
Płock, ul. Reja                      1
Radom, ul. Tochtermana               1
Siedlce, ul. Konarskiego             1
Skarżysko-Kamienna, MOBILNA          1
Starachowice, ul. Złota              1
Szczecin_Budzysza Wosia           1620
Toruń, ul. Dziewulskiego             1
Warszawa, al. Niepodległości         1
Warszawa, ul. Bajkowa                1
Warszawa, ul. Chrościckiego          1
Warszawa, ul. Kondratowicza          1
Warszawa, ul. Tołstoja               1
Warszawa, ul. Wokalna                1
Włocławek, ul. Okrzei                1
Zgierz, ul. Mielczarskiego           3
Łódź, ul. Czernika                   3
Łódź, ul. Gdańska                    3
Żyrardów, ul. Roosevelta             1
Name: PM25, dtype: int64

Dalej pozostały wartości brakujące które są jedyne w województwie - ich wartości uzupełniono poprzez średnią dla wszystkich stacji z tej samej data.

Pozostałe pojedyncze wartości uzupełniono poprzez interpolację.

Po tych zabiegach w zbiorze nie było już wartości brakujących.

data.loc[data['Nazwa stacji'] == 'Kalisz, ul. Wyszyńskiego', 'PM25'] = data['PM25'].fillna(data.groupby(data['Data'])['PM25'].transform('mean'))
data.loc[data['Nazwa stacji'] == 'Lublin, ul. Obywatelska', 'PM25'] = data['PM25'].fillna(data.groupby(data['Data'])['PM25'].transform('mean'))
data.loc[data['Nazwa stacji'] == 'Szczecin_Budzysza Wosia', 'PM25'] = data['PM25'].fillna(data.groupby(data['Data'])['PM25'].transform('mean'))
data['PM25'] = data['PM25'].interpolate()
brakujace_wartosci = data[data['PM25'].isnull()].groupby('Nazwa stacji')['PM25'].size()
brakujace_wartosci

Series([], Name: PM25, dtype: int64)

Ustawiono indeks danych jako kolumnę zawierającą datę i godzinę pomiaru.

data.set_index('Data', inplace=True)
data

	Nazwa stacji	PM25	Typ obszaru	Województwo	WGS84 φ N	WGS84 λ E
Data
2019-01-01 01:00:00	Duszniki-Zdrój	33.40530	miejski	DOLNOŚLĄSKIE	50.402645	16.393319
2019-01-01 02:00:00	Duszniki-Zdrój	13.80280	miejski	DOLNOŚLĄSKIE	50.402645	16.393319
2019-01-01 03:00:00	Duszniki-Zdrój	9.94056	miejski	DOLNOŚLĄSKIE	50.402645	16.393319
2019-01-01 04:00:00	Duszniki-Zdrój	6.75889	miejski	DOLNOŚLĄSKIE	50.402645	16.393319
2019-01-01 05:00:00	Duszniki-Zdrój	7.88722	miejski	DOLNOŚLĄSKIE	50.402645	16.393319
...	...	...	...	...	...	...
2019-12-31 20:00:00	Szczecin_Budzysza Wosia	15.59590	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294
2019-12-31 21:00:00	Szczecin_Budzysza Wosia	13.18250	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294
2019-12-31 22:00:00	Szczecin_Budzysza Wosia	14.37640	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294
2019-12-31 23:00:00	Szczecin_Budzysza Wosia	17.79410	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294
2020-01-01 00:00:00	Szczecin_Budzysza Wosia	20.86120	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294

499320 rows × 6 columns

Stworzono podzbiór zawierający dane zgrupowane dziennie (w celu lepszej czytelności wykresów). Wartość każdego dnia jest średnią wartością z danego dnia.

daily_avg_pm25 = data.groupby(['Nazwa stacji', pd.Grouper(freq='D')])['PM25'].mean().reset_index()
daily_avg_pm25

	Nazwa stacji	Data	PM25
0	Białystok, ul. Waszyngtona	2019-01-01	16.913043
1	Białystok, ul. Waszyngtona	2019-01-02	6.041667
2	Białystok, ul. Waszyngtona	2019-01-03	3.687500
3	Białystok, ul. Waszyngtona	2019-01-04	11.833333
4	Białystok, ul. Waszyngtona	2019-01-05	10.833333
...	...	...	...
20857	Żyrardów, ul. Roosevelta	2019-12-28	6.722500
20858	Żyrardów, ul. Roosevelta	2019-12-29	16.902917
20859	Żyrardów, ul. Roosevelta	2019-12-30	22.192083
20860	Żyrardów, ul. Roosevelta	2019-12-31	9.006250
20861	Żyrardów, ul. Roosevelta	2020-01-01	8.010000

20862 rows × 3 columns

station_groups = [daily_avg_pm25['Nazwa stacji'].unique()[i:i+5] for i in range(0, len(daily_avg_pm25['Nazwa stacji'].unique()), 5)]

for group in station_groups:
    group_data = daily_avg_pm25[daily_avg_pm25['Nazwa stacji'].isin(group)]
    plt.figure(figsize=(14, 10))
    for station in group:
        station_data = group_data[group_data['Nazwa stacji'] == station]
        plt.plot(station_data['Data'], station_data['PM25'], label=station)
    plt.title("PM2.5 dla różnych stacji")
    plt.xlabel("Data")
    plt.ylabel("PM2.5")
    plt.legend()
    plt.grid(True)

Stworozno wykresy liniowe ukazujące rozkład wartości PM2.5 w czasie dla wszystkich stacji. Widoczna jest pewna analogia, która ukazuje wysokie wahania i wartości stężeń w miesiącach chłodnych oraz niższe wartosci i wahania w ciepłym okresie.

srednie_pm25 = data.groupby('Województwo')['PM25'].mean()
srednie_pm25.plot(kind='bar', figsize=(10, 6), color='skyblue')
plt.title('Średnia wartość PM2.5 w zależności od województwa')
plt.xlabel('Województwo')
plt.ylabel('Średnia PM2.5')
plt.xticks(rotation=45)

(array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14]),
 [Text(0, 0, 'DOLNOŚLĄSKIE'),
  Text(1, 0, 'KUJAWSKO-POMORSKIE'),
  Text(2, 0, 'LUBELSKIE'),
  Text(3, 0, 'LUBUSKIE'),
  Text(4, 0, 'MAZOWIECKIE'),
  Text(5, 0, 'MAŁOPOLSKIE'),
  Text(6, 0, 'OPOLSKIE'),
  Text(7, 0, 'PODKARPACKIE'),
  Text(8, 0, 'PODLASKIE'),
  Text(9, 0, 'WARMIŃSKO-MAZURSKIE'),
  Text(10, 0, 'WIELKOPOLSKIE'),
  Text(11, 0, 'ZACHODNIOPOMORSKIE'),
  Text(12, 0, 'ŁÓDZKIE'),
  Text(13, 0, 'ŚLĄSKIE'),
  Text(14, 0, 'ŚWIĘTOKRZYSKIE')])

Średnie wartości dzienne w zależnosci od województwa są dosyć podobne. Wyrózniają się województwa podleskie oraz warmińsko-mazurskie (wartości najniższe) oraz mazowieckie (najwyższa wartość). Pozostałe rozkłądają się w przedziale średniej wartości 15-22.

top_wojewodztwa = data['Województwo'].value_counts().nlargest(5).index.tolist()
dane_top = data[data['Województwo'].isin(top_wojewodztwa)]

fig, axes = plt.subplots(nrows=len(top_wojewodztwa), ncols=1, figsize=(10, 6 * len(top_wojewodztwa)))

for i, wojewodztwo in enumerate(top_wojewodztwa):
    dane_wojewodztwo = dane_top[dane_top['Województwo'] == wojewodztwo]
    srednie_pm25 = dane_wojewodztwo.groupby(dane_wojewodztwo.index.hour)['PM25'].mean()
    srednie_pm25.plot(ax=axes[i], kind='line', marker='o', linestyle='-')
    axes[i].set_title(f'Średnia wartość PM2.5 dla {wojewodztwo}')
    axes[i].set_xlabel('Godzina')
    axes[i].set_ylabel('Średnia PM2.5')
    axes[i].set_xticks(range(24))
    axes[i].grid(True)

plt.tight_layout()

Następnie sprawdzono rozkład wartości w ciągu dnia dla pięciu najczęściej badanych województw.

data = data.interpolate(method='linear').dropna(axis = 1)
data.isna().sum()
fig, ax = plt.subplots(figsize=(12, 6))
acf = sm.tsa.acf(data['PM25'], nlags = 50)
plt.stem(range(len(acf)), acf)
plt.title(f'ACF dla PM2.5')
plt.xlabel('Opóźnienie')
plt.ylabel('Autokorelacja')
plt.grid(True)

C:\Users\Kowalski Krzysztof\AppData\Local\Temp\ipykernel_13264\1185333966.py:1: FutureWarning: DataFrame.interpolate with object dtype is deprecated and will raise in a future version. Call obj.infer_objects(copy=False) before interpolating instead.
  data = data.interpolate(method='linear').dropna(axis = 1)

Bazując na wykresie autokorelacji widoczna jest cykliczność skoków podobieństwa co dzień (24 godziny).

Tworzenie kolumn związanych z data

Dodano 3 kolumny związane z datą: dzień tygodnia (poniedziałek = 0, niedziela = 6), miesiąć oraz rok.

data['Day_of_week'] = data.index.dayofweek  
data['Month'] = data.index.month
data['Year'] = data.index.year
data

	Nazwa stacji	PM25	Typ obszaru	Województwo	WGS84 φ N	WGS84 λ E	Day_of_week	Month	Year
Data
2019-01-01 01:00:00	Duszniki-Zdrój	33.40530	miejski	DOLNOŚLĄSKIE	50.402645	16.393319	1	1	2019
2019-01-01 02:00:00	Duszniki-Zdrój	13.80280	miejski	DOLNOŚLĄSKIE	50.402645	16.393319	1	1	2019
2019-01-01 03:00:00	Duszniki-Zdrój	9.94056	miejski	DOLNOŚLĄSKIE	50.402645	16.393319	1	1	2019
2019-01-01 04:00:00	Duszniki-Zdrój	6.75889	miejski	DOLNOŚLĄSKIE	50.402645	16.393319	1	1	2019
2019-01-01 05:00:00	Duszniki-Zdrój	7.88722	miejski	DOLNOŚLĄSKIE	50.402645	16.393319	1	1	2019
...	...	...	...	...	...	...	...	...	...
2019-12-31 20:00:00	Szczecin_Budzysza Wosia	15.59590	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294	1	12	2019
2019-12-31 21:00:00	Szczecin_Budzysza Wosia	13.18250	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294	1	12	2019
2019-12-31 22:00:00	Szczecin_Budzysza Wosia	14.37640	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294	1	12	2019
2019-12-31 23:00:00	Szczecin_Budzysza Wosia	17.79410	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294	1	12	2019
2020-01-01 00:00:00	Szczecin_Budzysza Wosia	20.86120	miejski	ZACHODNIOPOMORSKIE	53.446761	14.507294	2	1	2020

499320 rows × 9 columns

Tworzenie kolumn bazując na średniej kroczącej i opoznieniach

data['1h_SR_K'] = data.groupby('Nazwa stacji')['PM25'].rolling(window=1, min_periods=1).mean().reset_index(level=0, drop=True)
data['1d_SR_K'] = data.groupby('Nazwa stacji')['PM25'].rolling(window=24, min_periods=1).mean().reset_index(level=0, drop=True)
data['3d_SR_K'] = data.groupby('Nazwa stacji')['PM25'].rolling(window=24 * 3, min_periods=1).mean().reset_index(level=0, drop=True)
data['7d_SR_K'] = data.groupby('Nazwa stacji')['PM25'].rolling(window=24 * 7, min_periods=1).mean().reset_index(level=0, drop=True)
data['30d_SR_K'] = data.groupby('Nazwa stacji')['PM25'].rolling(window=24 * 30, min_periods=1).mean().reset_index(level=0, drop=True)

data_fragment = data[(data.index.month == 12)]

fig, axes = plt.subplots(5, 1, figsize = [16, 12])
sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['PM25'], color = 'red', label = 'PM2.5', ax = axes[0])
sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['1h_SR_K'], color = 'blue', label = 'Średnia krocząca 1h', ax = axes[0])

sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['PM25'], color = 'red', label = 'PM2.5', ax = axes[1])
sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['1d_SR_K'], color = 'blue', label = 'Średnia krocząca 1d', ax = axes[1])

sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['PM25'], color = 'red', label = 'PM2.5', ax = axes[2])
sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['3d_SR_K'], color = 'blue', label = 'Średnia krocząca 3d', ax = axes[2])

sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['PM25'], color = 'red', label = 'PM2.5', ax = axes[3])
sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['7d_SR_K'], color = 'blue', label = 'Średnia krocząca 7d', ax = axes[3])

sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['PM25'], color = 'red', label = 'PM2.5', ax = axes[4])
sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['30d_SR_K'], color = 'blue', label = 'Średnia krocząca 30d', ax = axes[4])

<Axes: xlabel='Data', ylabel='PM25'>

data['1h_OP'] = data.groupby('Nazwa stacji')['PM25'].shift(1)
data['1d_OP'] = data.groupby('Nazwa stacji')['PM25'].shift(24)
data['3d_OP'] = data.groupby('Nazwa stacji')['PM25'].shift(24 * 3)
data['7d_OP'] = data.groupby('Nazwa stacji')['PM25'].shift(24 * 7)
data['30d_OP'] = data.groupby('Nazwa stacji')['PM25'].shift(24 * 30)

data_fragment = data[(data.index.month == 12)]

fig, axes = plt.subplots(5, 1, figsize = [16, 12])
sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['PM25'], color = 'red', label = 'PM2.5', ax = axes[0])
sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['1h_OP'], color = 'blue', label = 'Opoznienie 1h', ax = axes[0])

sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['PM25'], color = 'red', label = 'PM2.5', ax = axes[1])
sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['1d_OP'], color = 'blue', label = 'Opoznienie 1d', ax = axes[1])

sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['PM25'], color = 'red', label = 'PM2.5', ax = axes[2])
sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['3d_OP'], color = 'blue', label = 'Opoznienie 3d', ax = axes[2])

sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['PM25'], color = 'red', label = 'PM2.5', ax = axes[3])
sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['7d_OP'], color = 'blue', label = 'Opoznienie 7d', ax = axes[3])

sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['PM25'], color = 'red', label = 'PM2.5', ax = axes[4])
sns.lineplot(data_fragment, x = data_fragment.index, y = data_fragment['30d_OP'], color = 'blue', label = 'Opoznienie 30d', ax = axes[4])

<Axes: xlabel='Data', ylabel='PM25'>

2. Podział zbioru danych (5 punktów)

Podziel zbiór danych szeregów czasowych na zestawy treningowy i testowe, zachowując chronologię danych.

Pamiętaj, że posiadasz dane o charakterystyce czasowo-przestrzennej.

Wykorzystaj ten podział do wstępnej oceny działania modelu.

W celu podziału danych u względnieniem czasu przed podziałem posortowano dane oraz dokonano kategoryzacji danych innych niż liczbowe.

data.sort_index(inplace = True)

X = data.drop(columns=['PM25', '1h_SR_K'])  #'1h_OP'
y = data['PM25']  

encode_cols = ['Nazwa stacji', 'Typ obszaru', 'Województwo']

for col in encode_cols:
    labelencoder = LabelEncoder()
    X[col] = labelencoder.fit_transform(X[col])
X

	Nazwa stacji	Typ obszaru	Województwo	WGS84 φ N	WGS84 λ E	Day_of_week	Month	Year	1d_SR_K	3d_SR_K	7d_SR_K	30d_SR_K	1h_OP	1d_OP	3d_OP	7d_OP	30d_OP
Data
2019-01-01 01:00:00	4	0	0	50.402645	16.393319	1	1	2019	13.000000	13.000000	13.000000	13.000000	NaN	NaN	NaN	NaN	NaN
2019-01-01 01:00:00	50	0	12	51.856692	19.421231	1	1	2019	76.730600	76.730600	76.730600	76.730600	NaN	NaN	NaN	NaN	NaN
2019-01-01 01:00:00	32	0	7	50.040675	22.004656	1	1	2019	42.034031	42.034031	42.034031	42.034031	NaN	NaN	NaN	NaN	NaN
2019-01-01 01:00:00	28	0	7	49.784339	22.756239	1	1	2019	38.781605	38.781605	38.781605	38.781605	NaN	NaN	NaN	NaN	NaN
2019-01-01 01:00:00	48	0	3	51.799722	16.317500	1	1	2019	64.075200	64.075200	64.075200	64.075200	NaN	NaN	NaN	NaN	NaN
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
2020-01-01 00:00:00	31	1	7	49.546539	21.851006	2	1	2020	9.960000	17.129028	15.910893	22.857958	13.6000	22.00000	13.60000	6.1000	17.4000
2020-01-01 00:00:00	53	0	12	51.758050	19.529786	2	1	2020	11.926698	29.433478	17.405045	25.337234	9.0000	21.00000	9.00000	44.0000	37.0000
2020-01-01 00:00:00	32	0	7	50.040675	22.004656	2	1	2020	9.446264	14.457141	13.033345	19.184583	15.2866	18.55700	29.71530	17.7406	15.9014
2020-01-01 00:00:00	13	0	14	50.878998	20.633692	2	1	2020	13.712149	17.139649	13.150660	23.136434	8.5700	25.56000	5.44000	20.7500	55.7900
2020-01-01 00:00:00	36	0	11	53.446761	14.507294	2	1	2020	8.566667	16.024861	15.851786	20.419736	17.7941	7.33104	9.37627	16.9329	18.6052

499320 rows × 17 columns

Model z opóźnieniem 1h i średnią krocząca 1h - taki model osiągnie lepsze wyniki natomiast będzię nadawać się do przewidywania najbliższej przyszłości. Dane podzielono w sposób: 20% dane testowe oraz 80% dane treningowe.

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.20, shuffle = False)

fig, axes = plt.subplots(1,1, figsize = [12, 8])
plt.plot(y_train, color = 'red')
plt.plot(y_test, color = 'blue')

[<matplotlib.lines.Line2D at 0x158b102bfd0>]

Część 2: Tworzenie modelu do predykcji (20 punktów)

1. Inicjalizacja i trenowanie modelu (10 punktów)

W tym ćwiczeniu do predykcji będziemy używaż algorytmu XGB (Extreme Gradient Boosting).

Zainicjuj i wytrenuj model na danych uczących.

Napisz plusy i minusy wykorzystania tego algorytmu w kontekście predykcji szeregów czasowych w formie zwięzłej tabeli

XGB Boost plusy:

• posiada wysoką dokładność predykcji dzięki wykorzystywaniu skomplikowanych zależności nieliniowych,
• dane związane z szeregami czasowymi, często posiadają wartości brakujące, XGB Boost posiada wbudowane mechanizmy radzenia sobie z takimi wartościami,
• algorytm posiada techniki L1 i L2 (regularyzacji), które pozalajaą zapobiec przetrenowaniu się modelu (częsty bład w przypadku szeregów czasowych).

XGB Boost minusy:

• proces przygotowania danych i konfiguracji modelu może być czasochłonny i problematyczny,
• algorytm posiada wiele hiperparametrów, których wartość należy dostroić, aby zwiększych efekty,
• XGB Boost może być zasobobieżny, za równo pod kontem pamięci jak i złożoności obliczeniowej.

model = xgb.XGBRegressor()
model.fit(X_train, y_train)
y_pred = model.predict(X_test)

mse = mean_squared_error(y_test, y_pred)
print("Błąd średniokwadratowy (MSE):", mse)

r2 = r2_score(y_test, y_pred)
print("Współczynnik determinacji (R-squared):", r2)

Błąd średniokwadratowy (MSE): 77.01389925010608
Współczynnik determinacji (R-squared): 0.770628732017302

2. Tuning hiperparametrów (10 punktów)

Wykonaj tuning hiperparametrów poprzez trening i walidację krzyżową w celu znalezienia optymalnych wartości hiperparametrów.

Napisz czym są hiperparametry i dlaczego wymagają tuningu.

Hiperparametry to parametry algorytmu, które są ustawiane przed rozpoczęciem procesu uczenia maszynowego i które nie są optymalizowane podczas treningu modelu. Są to parametry, które kontrolują zachowanie algorytmu i strukturę modelu. Wpływają na jego wydajność i skuteczność. Tuning ich pozwala na lepsze dobranie ich co usprawni proces uczenia, a dzięki temu wynik całego działania algorytmu. Ponieważ algorytm XGB Boost posiada wiele hiperparametrów ich dostosowanie pozwoli osiągnąć lepsze wyniki dla badanego problemu.

param_grid = {
    'max_depth': [3, 5, 7],
    'learning_rate': [0.1, 0.01, 0.001],
    'n_estimators': [100, 200, 300]
}
grid_search = GridSearchCV(estimator=xgb.XGBRegressor(), param_grid=param_grid, cv=3, scoring='neg_mean_squared_error')
grid_search.fit(X_train, y_train)

print("Najlepsze parametry:", grid_search.best_params_)
print("Najlepszy wynik:", grid_search.best_score_)

Najlepsze parametry: {'learning_rate': 0.1, 'max_depth': 3, 'n_estimators': 100}
Najlepszy wynik: -69.81401232093168

best_model = xgb.XGBRegressor(**grid_search.best_params_)
best_model.fit(X_train, y_train)
y_pred = best_model.predict(X_test)

mse = mean_squared_error(y_test, y_pred)
print("Błąd średniokwadratowy (MSE):", mse)

r2 = r2_score(y_test, y_pred)
print("Współczynnik determinacji (R-squared):", r2)

Błąd średniokwadratowy (MSE): 62.87115058269354
Współczynnik determinacji (R-squared): 0.812750222114436

index = range(len(y_test))
plt.figure(figsize=(10, 6))
plt.plot(index, y_test, label='Rzeczywiste wartości', color='red')
plt.plot(index, y_pred, label='Przewidziane wartości', color='blue')
plt.xlabel('Indeks danych testowych')
plt.ylabel('Wartości PM25')
plt.title('Porównanie rzeczywistych i przewidywanych wartości')
plt.legend()
plt.grid(True)

Część 3: Backtesting (20 punktów)

1. Podział na okna czasowe (10 punktów)

Podziel dane na wiele okien czasowych, na przykład poprzez użycie walk-forward lub roll-forward validation lub gap-roll-forward validation, aby uzyskać rzeczywistą ocenę wydajności modelu na przyszłych danych.

Napisz na czym polega ten sposób walidacji i po co się go wykonuje.

Walidacja na wielke okien czasowych jest kluczową techniką w modelowaniu szeregów czasowych. Zapewnia ona realistyczną ocenę modelu w kontekście zmieniających się danych. Polega ona na iteracyjnym trenowaniu i testowaniu modelu na różnych podzbiorach danych, przesuwając okno czasowe z każdą iteracją. Dzięki tej metodzie możemy lepiej zrozumieć, jak model radzi sobie w różnych warunkach czasowych i uniknąc błędów które mogą być powodowane przez tradycyje metody walidacji.

tscv = TimeSeriesSplit(n_splits = 8)

2. Ocena wydajności (10 punktów)

Przeprowadź backtesting modelu, czyli ocenę wydajności modelu na danych historycznych, które nie zostały użyte do trenowania.

Czym rózni się backtesting od testowania wydajności modelu dla danych niebędących szeregami czasowymi?

Główne różnice między backtestingiem, a testowaniem wydajności modelu dla danych niebędących szeragami czasowymi to?

• w backtestingu kolejność danych jest kluczowa i musi być zachowana, aby można było symulować przyszłe prognozy. W klasycznym kolejność danych nie ma znaczenia, a podziały mogą być losowe,
• backtesting odzwierdziedla rzeczywiste warunki pracy modelu w przyszłości, co lepiej sprawdza się w analizie szeregów czasowych,
• backtesting korzysta z przesuwających się okien czasowych (walk-forward, roll-forward lub gap-roll-forward validation), gdy klasycznie zazwyczaj korzysta się z metod losowego podziału takich jak crosswalidacja,
• dane czasowe possiadają naturalne zależności między kolejnymi obserwacjami, co musi być uwzgłędnione w procesie walidacji. W danych klasycznach obserwację są zazwyczaj niezależne.

mse_scores = []
r2_scores = []

model = xgb.XGBRegressor(**grid_search.best_params_)
for train_index, test_index in tscv.split(X_train):
    X_train_fold, X_test_fold = X_train.iloc[train_index], X_train.iloc[test_index]
    y_train_fold, y_test_fold = y_train.iloc[train_index], y_train.iloc[test_index]
    
    model.fit(X_train_fold, y_train_fold)
    
    y_pred = model.predict(X_test_fold)
    
    mse = mean_squared_error(y_test_fold, y_pred)
    mse_scores.append(mse)
    
    r2 = r2_score(y_test_fold, y_pred)
    r2_scores.append(r2)

mean_mse = np.mean(mse_scores)
mean_r2 = np.mean(r2_scores)

print("Mean MSE:", mean_mse)
print("Mean R^2:", mean_r2)

Mean MSE: 29.682145150692634
Mean R^2: 0.8011858514966412

Część 4: Analiza istotności cech (20 punktów)

1. Wyliczenie istotności cech (5 punktów)

Wykorzystaj wbudowane narzędzia do wyliczenia istotności cech w modelu (Feature Importance).

feature_importance = model.feature_importances_
columns = list(X_train.columns)
for i, importance in enumerate(feature_importance):
    print(f"Istotność cechy {columns[i]}: {importance}")

Istotność cechy Nazwa stacji: 0.03198323771357536
Istotność cechy Typ obszaru: 0.00031548592960461974
Istotność cechy Województwo: 0.007906577549874783
Istotność cechy WGS84 φ N: 0.002865417627617717
Istotność cechy WGS84 λ E: 0.0017405450344085693
Istotność cechy Day_of_week: 0.0017863311804831028
Istotność cechy Month: 0.016754381358623505
Istotność cechy Year: 0.0
Istotność cechy 1d_SR_K: 0.00785788893699646
Istotność cechy 3d_SR_K: 0.02538600005209446
Istotność cechy 7d_SR_K: 0.006205314304679632
Istotność cechy 30d_SR_K: 0.0019663472194224596
Istotność cechy 1h_OP: 0.853451132774353
Istotność cechy 1d_OP: 0.019817249849438667
Istotność cechy 3d_OP: 0.008371957577764988
Istotność cechy 7d_OP: 0.004258961882442236
Istotność cechy 30d_OP: 0.009333128109574318

plt.figure(figsize=(10, 6))
plt.bar(columns, feature_importance)
plt.xlabel('Numer cechy')
plt.ylabel('Istotność')
plt.title('Istotność cech')
plt.xticks(columns, rotation = 90)
plt.grid(axis='y')

2. Interpretacja wyników (15 punktów)

Zinterpretuj istotność cech w kontekście dostępnych informacji - artykułów naukowych, dokumentacji biblioteki itp.

W powyższym modelu posiadającym kolumne bazującą na opóznieniu wysokości 1h widoczny jest ogromny wpływ tej kolumny na działanie modelu. Bazując na otrzymanych wartościach metryk testowych model ten bardzo dobrze przewiduje wartości PM2.5. W celu próby stworzenia bardziej uniwersalnego modelu usunięto cechę ukazującą opóźnienie jednogodzinne.

data.sort_index(inplace = True)

X = data.drop(columns=['PM25', '1h_SR_K',  '1h_OP'])
y = data['PM25']  

encode_cols = ['Nazwa stacji', 'Typ obszaru', 'Województwo']

for col in encode_cols:
    labelencoder = LabelEncoder()
    X[col] = labelencoder.fit_transform(X[col])
X

	Nazwa stacji	Typ obszaru	Województwo	WGS84 φ N	WGS84 λ E	Day_of_week	Month	Year	1d_SR_K	3d_SR_K	7d_SR_K	30d_SR_K	1d_OP	3d_OP	7d_OP	30d_OP
Data
2019-01-01 01:00:00	4	0	0	50.402645	16.393319	1	1	2019	13.000000	13.000000	13.000000	13.000000	NaN	NaN	NaN	NaN
2019-01-01 01:00:00	50	0	12	51.856692	19.421231	1	1	2019	76.730600	76.730600	76.730600	76.730600	NaN	NaN	NaN	NaN
2019-01-01 01:00:00	32	0	7	50.040675	22.004656	1	1	2019	42.034031	42.034031	42.034031	42.034031	NaN	NaN	NaN	NaN
2019-01-01 01:00:00	28	0	7	49.784339	22.756239	1	1	2019	38.781605	38.781605	38.781605	38.781605	NaN	NaN	NaN	NaN
2019-01-01 01:00:00	48	0	3	51.799722	16.317500	1	1	2019	64.075200	64.075200	64.075200	64.075200	NaN	NaN	NaN	NaN
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
2020-01-01 00:00:00	31	1	7	49.546539	21.851006	2	1	2020	9.960000	17.129028	15.910893	22.857958	22.00000	13.60000	6.1000	17.4000
2020-01-01 00:00:00	53	0	12	51.758050	19.529786	2	1	2020	11.926698	29.433478	17.405045	25.337234	21.00000	9.00000	44.0000	37.0000
2020-01-01 00:00:00	32	0	7	50.040675	22.004656	2	1	2020	9.446264	14.457141	13.033345	19.184583	18.55700	29.71530	17.7406	15.9014
2020-01-01 00:00:00	13	0	14	50.878998	20.633692	2	1	2020	13.712149	17.139649	13.150660	23.136434	25.56000	5.44000	20.7500	55.7900
2020-01-01 00:00:00	36	0	11	53.446761	14.507294	2	1	2020	8.566667	16.024861	15.851786	20.419736	7.33104	9.37627	16.9329	18.6052

499320 rows × 16 columns

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.20, shuffle = False)

fig, axes = plt.subplots(1,1, figsize = [12, 8])
plt.plot(y_train, color = 'red')
plt.plot(y_test, color = 'blue')

[<matplotlib.lines.Line2D at 0x158b5f88cd0>]

model = xgb.XGBRegressor()
model.fit(X_train, y_train)
y_pred = model.predict(X_test)

mse = mean_squared_error(y_test, y_pred)
print("Błąd średniokwadratowy (MSE):", mse)

r2 = r2_score(y_test, y_pred)
print("Współczynnik determinacji (R-squared):", r2)

Błąd średniokwadratowy (MSE): 311.11213723639395
Współczynnik determinacji (R-squared): 0.07341160365126487

Bez tej cechy model posiada znacznie gorsze statystyki.

param_grid = {
    'max_depth': [3, 5, 7],
    'learning_rate': [0.1, 0.01, 0.001],
    'n_estimators': [100, 200, 300]
}
grid_search = GridSearchCV(estimator=xgb.XGBRegressor(), param_grid=param_grid, cv=3, scoring='neg_mean_squared_error')
grid_search.fit(X_train, y_train)

print("Najlepsze parametry:", grid_search.best_params_)
print("Najlepszy wynik:", grid_search.best_score_)

Najlepsze parametry: {'learning_rate': 0.1, 'max_depth': 5, 'n_estimators': 100}
Najlepszy wynik: -217.39854364121285

best_model = xgb.XGBRegressor(**grid_search.best_params_)
best_model.fit(X_train, y_train)
y_pred = best_model.predict(X_test)

mse = mean_squared_error(y_test, y_pred)
print("Błąd średniokwadratowy (MSE):", mse)

r2 = r2_score(y_test, y_pred)
print("Współczynnik determinacji (R-squared):", r2)

Błąd średniokwadratowy (MSE): 282.02101492499406
Współczynnik determinacji (R-squared): 0.16005398478737332

index = range(len(y_test))
plt.figure(figsize=(10, 6))
plt.plot(index, y_test, label='Rzeczywiste wartości', color='red')
plt.plot(index, y_pred, label='Przewidziane wartości', color='blue')
plt.xlabel('Indeks danych testowych')
plt.ylabel('Wartości PM25')
plt.title('Porównanie rzeczywistych i przewidywanych wartości')
plt.legend()
plt.grid(True)

tscv = TimeSeriesSplit(n_splits = 8)

mse_scores = []
r2_scores = []

model = xgb.XGBRegressor(**grid_search.best_params_)
for train_index, test_index in tscv.split(X_train):
    X_train_fold, X_test_fold = X_train.iloc[train_index], X_train.iloc[test_index]
    y_train_fold, y_test_fold = y_train.iloc[train_index], y_train.iloc[test_index]
    
    model.fit(X_train_fold, y_train_fold)
    
    y_pred = model.predict(X_test_fold)
    
    mse = mean_squared_error(y_test_fold, y_pred)
    mse_scores.append(mse)
    
    r2 = r2_score(y_test_fold, y_pred)
    r2_scores.append(r2)

mean_mse = np.mean(mse_scores)
mean_r2 = np.mean(r2_scores)

print("Mean MSE:", mean_mse)
print("Mean R^2:", mean_r2)

Mean MSE: 151.57003149244986
Mean R^2: 0.20378762897091302

feature_importance = model.feature_importances_
columns = list(X_train.columns)
for i, importance in enumerate(feature_importance):
    print(f"Istotność cechy {columns[i]}: {importance}")

Istotność cechy Nazwa stacji: 0.04847313091158867
Istotność cechy Typ obszaru: 0.045297861099243164
Istotność cechy Województwo: 0.015866482630372047
Istotność cechy WGS84 φ N: 0.05105580762028694
Istotność cechy WGS84 λ E: 0.027293642982840538
Istotność cechy Day_of_week: 0.03461988642811775
Istotność cechy Month: 0.06915411353111267
Istotność cechy Year: 0.0
Istotność cechy 1d_SR_K: 0.23705363273620605
Istotność cechy 3d_SR_K: 0.013916288502514362
Istotność cechy 7d_SR_K: 0.03167867660522461
Istotność cechy 30d_SR_K: 0.019020305946469307
Istotność cechy 1d_OP: 0.29580432176589966
Istotność cechy 3d_OP: 0.039913661777973175
Istotność cechy 7d_OP: 0.03651047870516777
Istotność cechy 30d_OP: 0.03434166684746742

plt.figure(figsize=(10, 6))
plt.bar(columns, feature_importance)
plt.xlabel('Numer cechy')
plt.ylabel('Istotność')
plt.title('Istotność cech')
plt.xticks(columns, rotation = 90)
plt.grid(axis='y')

WNIOSKI

• W danych widoczna była dzienna zzmienność danych - wartości zanieczyszczenia PM2.5 w ciągu dnia są wysokie w godzinach nocnych, dalej często występujące w Polsce piece (np. typu koza) wymagają wieczornego dorzucenia, aby utrzymać ciepło co widoczne jest w wykresie - wzrost wartości w godzinach, gdy ludzie zazwyczaj kładą się spać oraz delikatny pik po wstaniu (godzina 6-7). Podczas gdy większość osób jest wp pracy wartość stężenia jest niższa.
• Widoczny jest rónież wpływ temperatury powietrza (miesiące zimne, a ciepłe) na wielkość zanieczyszenia oraz wahania tej wielkości, co również może mieć związek z ogrzewaniem.
• Algorytm XGB Boost pozwala osiągać ciekawe i dobre wyniki działania na danych z czynnikiem czasowym.
• Tuning hiperparametrów pozwala dobrać je tak, aby końcowy model był jeszcze lepszy niż otrzymany.
• Backtesting jest ważny podczas walidacji modelu szeregu czasowego. Pozwala na przetestowanie go w różnych warunkach podobnych do rzeczywistych.
• Sprawdzenie istotności cech uwypukla ważność każdej z cech, co może pomóc wykryć błędne założenia lub wyciągnąc dodatkowe wnioski z danych.
• Stworzenie dwóch różnych modeli (jeden posiadający cechę z opóznieniem godzinowym) pozwoliło mi zwrócic uwagę jak dane te (o czym świadczył wykres autokorelacji, a po wykonaniu modelu wykres istotności cech) znacznie wpływają na model. Drugi z modeli bez tej cechy posiada wyrażnie gorsze zdolności predykcyjne i jest po prostu gorszym modelem, natomiast może on zostać wykorzystany do predykcji dalszych wartości w przyszłosci niż pierwszy ze stworzonych modeli. Co po raz kolejny pokazuje jak dobór i tworzenie cech ma znaczoncy wpływ na końcowy wygląd modelu.
• Dane, które posiadają informacje o stężeniu zanieczyszczenia PM2.5 dobrze nadają się do testowania na nich metod predykcji szeregów czasowych oraz pozwalają osiągnąć dobre wyniki analizy.

Część ogólna (10 punktów)

1. Zgodność kodu z PEP-8 oraz elegancja kodu

2. Ogólna interpretacja wyników

3. Skuteczność i precyzja

Oceniając wyniki backtestingu, skup się na skuteczności i precyzji predykcji modelu.

4. Kreatywność i innowacyjność

Oceniając analizę istotności cech oraz tworząć cechy poszukaj kreatywnych i innowacyjnych sposobów działania.

5. Strona edytorska

Do oceny przekaż tylko finalny kod oraz interpretację. Pamiętaj o staranności edytorskiej.

Oddanie projektu w terminie określonym w systemie (10 punktów)